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ABSTRACT 

We measure the monopole moment of the three-point correlation function on scales 
lh- x Mpc - 70h- x Mpc in the Two degree Field Galaxy Redshift Survey (2dFGRS). 
Volume limited samples are constructed using a series of integral magnitudes bins 
between M = —18 ... — 22. Our measurements with a novel edge-corrected estimator 
represent most, if not all, three-point level monopole or angular averaged information 
in the catalogue. We fit a perturbative non-linear bias model to a joint data vector 
formed from the estimated two- and three-point correlation functions. Two different 
models are used: an analytic model based on Eulerian perturbation theory including 
bias and redshift distortions, and a phenomenological bias model based on the direct 
redshift space measurements in the large Virgo simulations. To interpret the clustering 
results, we perform a three parameter Gaussian maximum likelihood analysis. In the 
canonical —21 ~ —20 volume limited sample we find erg = 0.93j^Q "jj , b = 1.04^qq|, 

and 62 = —0.06^0 001- Our estimate of eg, is robust across the different volume limited 
samples constructed. These results, based solely on the large scale clustering of galax- 
ies, are in excellent agreement with previous analyses using the Wilkinson Anisotropy 
Probe: this is a spectacular success of the concordance model. We also present two- 
parameter fits for the bias parameters, which are in excellent agreement with previous 
findings of the bias evolution in the 2dFGRS. 
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1 INTRODUCTION 

Statistical analyses of t he Two Degree Fie ld Galaxy Red- 
shift Survey (2dFGRS) (IColfess et al]l200lft have propelled 
significant progress in high precision cosmology. For in- 
stance, measurements of the two-point correlation function 
and the power spectrum have provided tight constraints on 
the theories of struc t ure formation (e.fi. |Percival^^alJ[200ll: 
iNorberg et aill200ll: lHawkins et alll2003t ICole et alJ^OOsf . 
The large volume and high quality of the 2dFGRS encourage 
further studies of higher order statistics. Such investigations 
provide information on the Gaussianity of the small initial 
dark matter fluctuations, the emergence of non- Gaussianity 
through non-linear gravitational effects, and even the murky 
physical processes of ga laxy formati on. The latter might 
manifest itself as "bias" (Kaiser 1984), where the clustering 
of galaxies might be statistically different from that of the 
dark matter. Higher than second order statistics provide the 
only tool with which to separate these effects from the grav- 
itational amplification and initial conditions. The resulting 
constraints on the bias are interesting in their own right, 



and they provide the new avenues to ultimately constrain 
cosmological parameters. 

Third order statistics represent the first non-trivial step 
in the perturbative understanding of non-Gaussianity. In- 
deed, to date numerous works have been devoted mainly 
to the measurement an d understanding o f the third or- 
der statistics of 2dFGRS. IVerde et al.l i2002F) estimated bias 
parameters from their bispectrum measurement at wave 
lengths 0.1 < k < 0.5h/Mpc. Three point correlation func- 
tion £ of an early releas e d 2dF sample (2dF100k) is mea- 
sured by Ijing fc Bornerl <| 2004l) focusing on empirical for- 
mula fit for c". IWane^tHuT(l27)04h measured three point cor- 
relation function for the 2dFGRS on small scales to test their 
conditional luminosity function model jointly with the halo 
model. As an alter native to the three-p oint correlation func- 
tion or bispectrum. ICroton et all (l2004h calculated moments 
of counts in cells, or averaged iV-point correlation functions, 
on 1 — 9h~ Mpc scales, and estimated relative bias param- 
eters. 

A common thread in previous measurements was that 
they focused on either relatively small scales and/or a set 
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of hand picked subset of triangular configurations, which 
characterize three-point statistics. With the most natural 
parametrizations predominantly used in the past, such as 
the three sides of a triangle C( r ii r 2, r$), or £(ri,r2,#) with 
6 the angle between ri and f2, it is both computationally 
burdensome, and conceptually difficult to measure and in- 
terpret thj-«^point_statistics in a large dynamic range. Re- 
cently, ISzapudil <l2004h has shown that a multipole expan- 
sion motivated by rotational invariance helps substantially 
with this "combinatorial expansion" of parameters. It was 
demonstrated that it is most efficient to expand spherically 
symmetric functions of two unit vectors into L = bipo- 
lar spherical harmonics. In turn, the recipe boils down to 
expand £ into Legendre polynomials Pi(cos8) 



C(n,r- 2 ,( 



^9/4-1 



(i) 



ISzapudi J2004) has demonstrated that the first few multi- 
poles, often up to t — 2, concentrate most of the useful 
three-point level information. While one still needs to con- 
sider the scale dependence over r\ and r2, the configuration 
space effectively becomes two-dimensional. 

In addition to the conceptual simplification, the 
above parametrization sug gests new algor ithms to calcu- 
late three-point funct ions llSzapudil 1200 5a') . In particular, 
iPan fc S zapudi| (|2005]) adapted the fully edge corrected es- 
timator of ISzapudi fc Szalavl ( f 998) for the £ = monopole 



moment, and have demonstrated a simple and fast JV 2 algo- 
rithm for its calculation. Edge correction is a major advan- 
tage over other three-point statistics which can be calculated 
reasonably fast, most notably the bispectrum and moments 
of counts in cells. This means that the resulting estimates of 
the monopole of the three-point function are expected to be 
more robust against complicated geometry of the window, 
cut out holes etc, than other previously used measures. Since 
real surveys, such as the 2dFGRS have complicated spatial 
structure, edge effect correction is a must when approaching 
large scales. 

The monopole moment captures all information about 
the amplitude of the three-point correlation function; all 
other multipoles provide information on the shape. Besides 
the fact that the monopole moment is the lowest order, thus 
the simplest to measure and interpret in the multipole se- 
ries, it also has a simple transformation under bias, and we 
have a r elatively accurate und erstanding of its redshift dis- 
tortions iPan fc Szapudill2005l) . These properties single out 
the monopole moment as a principal candidate among the 
three-point statistics for practical applications. 

In this paper we set out to harvest the fruits of recent 
theoretical developments, and measure and interpret the 
monopole moment of the three-point function in the 2dF- 
GRS. The interpretation of three-point st atist i cs in terms 
of bias parameters was put forward by iFrvl £l994) . His 
method was later perfected to include maximum likelihood 



of ratio statistics jMatarrese et alJll997l: 


Verde et alJ 
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Scoccimarrol l2000l: iFeldman et al.l 1200 ft 


Verde et alJ 


2002 


Gaztanaga fc Scoccimarroll2005r). We develop a novel joint 



maximum likelihood technique using both two- and three- 
point statistics (as opposed to ratio statistics) for simultane- 
ous estimation of bias coefficients and cosmological param- 
eters, such as as- We estimate covariances in the data using 



mock surveys, and constrain the parameters of our theory 
in a Gaussian maximum likelihood context with scales up 
to 140h _1 Mpc entering into the analysis. Even though £0 
is only the first in the series of multipoles, we will see that 
it contains invaluable, hitherto untouched information on 
cosmology and bias. 

The next section outlines our method of estimating £0 
in volume limited subsamples of 2dFGRS; section 3 details 
the theoretical framework for the interpretation of the data 
in terms of bias and cosmological parameters; the resulting 
constraints are presented in section 4; discussion and sum- 
mary follows in section 5. 



2 MEASUREMENT OF Co 
2.1 The data set 

In order to estimate three-point correlation functions, we 
constructed volume limited samples from the 2dFGRS fi- 
nal data release spectroscopic catalogue (the 2dF230k, 
IColless et ail 12003ft with 221414 galaxies of good redshift 
quality Q > 3 l|Colless et alJl200lfr . Excluding the ancillary 
random fields leaves us two large contiguous volumes: one 
near the South Galactic Pole (SGP) covering approximately 
-37°.5 < 5 < -22°.5, 21 h 40 m < a < 3 h 40 m and the other 
one around the North Galactic Pole (NGP) defined roughly 
by -715 < S < 2°.5, 9 h 50 m < a < 14 h 50 m . The parent 
sample was further restricted by completeness / > 0.7, and 
apparent magnitudes limits in photometric bj band with 
bright cut of mtj = 15 and faint cut of median value of 
~ 19.3 with certain small variation as specified by masks 
jColless et al.ll200ll) . 

Volume limited subsamples are built from the parent 
sample by selecting galaxies in specified absolute magnitude 
ranges. These were calculated with k + e correction as in 
iNorberg et alJ (2002). The most important properties of the 
result SGP and NGP samples are summarized in Table 
For our measurements, the NGP and SGP were combined 
together to achieve the highest possible volume. 



2.2 Estimation of the correlation function 

Two point correlation functions are measured with the 
lLandv fc Szalavl lll993f) estimator 



DD - 2DR + RR 
RR 



(2) 



Here D stands for data and R for points selected from ran- 
dom catalogues. These were created according to the exact 
geometry and completeness masks of the subsamples, with 
20 times as many random points as the number of galaxies 
N B . 

The three point correlat ion function also has s imilar 
minimum variance estimator iSzapudi fc Szalavl ll998l 



C 



DDD - 3DDR + 3DRR - RRR 
RRR 



(3) 



This estimator has been shown bv lKavo et al.l i2004f) to be 
more accurate than other estimators currently in use. Be- 
cause of the volume limited samples, no additional weighting 
was necessary. 

For the present analyses, we measured the monopole 
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Table 1. Volume limited subsamples of 2dFGRS. Comoving distance d is calculated from redshifts z assuming ACDM universe with 
fl A = 0.7, tt m = 0.3. 



M bj - 


5 log 10 h 




Zmax 


drain 




N^gp/ n ngp 


n(10" 3 /i 3 Mpc~ 3 ) 


-18 - 


-17 


0.0131 


0.0575 


39.0 


170.0 


4046/3192 


12.97 


-19 - 


-18 


0.0205 


0.087 


61.2 


255.7 


11935/9625 


11.35 


-20 - 


-19 


0.0320 


0.129 


95.2 


374.9 


23595/17729 


6.922 


-21 - 


-20 


0.0495 


0.186 


146.6 


532.9 


18081/12499 


1.798 


-22 — 


-21 


0.0754 


0.261 


222.2 


735.7 


4095/2113 


0.140 



moment or angular ave rage of the three-point correlation 
function iSzapudill20ol) 



(4) 



Co(Vi, T2) = 27r / dcos 8( (ri, rz, 0) 



Since this is an angular averaged quantity, the estimator 
of Eq. |21 can be realized with setting up bins of shells 
around center points. These are essentially neighbour counts 
in shells, a nd can be realized with a simple iV 2 estimator put 
forward in lPan fc SzapudU i2005l) . Explicitly, if the number 
of galaxies around galaxy i in bin {r°,r-^) is m(ri), in bin 
{r l 2 °,rf) is ni(r 2 ), the DDD in Eq.EJreads 



S i j 1 "i('-i)("»('-2)-i) 

JV B (iV g -l)(JV g -2) 



if n =fc r2 
if n = r 2 



(5) 



This estimator has no Poisson noise bias as there is no 
overlap between the configurations. The lack of shot noise 
bias, and the precise edge correction are major techni- 
cal advantages over measuring the bispectrum via direc t 
Fourier transform (e.g. lScoccimarrol2000HVerde et all2002l). 
or moments of counts in cells (e.p ; .. [Szapudi fc Colombil 
Il996l: ISzapudilll99l iBernardeau et alJl200l and references 
therein). 

Redshift space measurements were performed in 29 log- 
arithmic bins between 0.7 — 70h _1 Mpc for the two-point cor- 
relation function, and all corresponding pairs for CofYi,?^), 
altogether 29 x 30/2 = 435 bins. All measurements were 
repeated in the available 22 mock 2dF surveys, as well as 
eight equal sub cubes selected from the Virgo VLS ACDM 
simulations iMacfarland et aT]ll998t) . 



3 INTERPRETATION 
3.1 Theoretical Framework 

To interpret the clustering present in the 2dF, we use two 
models: theoretical model for large scales, and a phenomeno- 
logical model for intermediate scales. 

The t heoretical model uses E ulerian perturbation the- 
ory (e.g.. IBernardeau et aljl2002l and references therein) 
to calculate the real space dark matter three-poi nt corre- 
lation f unction; specifically, we use the formulae m lSzapudil 
J2005bl) . 

Our bias model is motivate d by the usual perturbat ive 
expansion S g = b8 + &2<5/2 + . . . (IFrv fe Gaztanagalll993i) 
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Figure 1. We display our measurements of the two- and three- 
point correlation functions £ (red solid circles) and (o (blue solid 
squares) in the canonical M=— 21 ~ —20 volume limited sam- 
ple. Error bars have been calculated from 22 2dF mocks. Solid 
line is our phenomenological model with parameters in Table 151 
dashed lines are theoretical model with parameters in Table |^] 
Inset plotted on a linear scale for on large scales. 



where 6 and 62 are the linear and non-linear bias factors, 
(...)o denotes angular averagi ng, and fa, fa are redshift dis- 
tortion enhancement factors (Ka iserl Il987l : lHamilto'nlll99ftf) 



_ 5(2520 + 3360/ + 1260/ 2 + 9/ 3 - 14/ 4 ) 7 2 

» ^ no /1 rr i -i r\ £ i o £0 \ 9 ^ 



(7) 



where / 



98(15 + 10/ + 3/ 2 ) 2 
Q°.6/b. The third order f^ is obtained by 



IPan fc Sza pudi ( 2005) th rough angle avera gi ng perturbation 
theory results of IScoccimarro et alJ jl999^|Pan fc Szapudil 
( 2005) have shown this this is a good approximation to angle 
averaged three-point quantities with deviation at 5% level. 

Our phenomenological model is obtained by using di- 
rectly our measured Virgo VLS redshift space two and three- 
point correlation functions in Equation |S| (with fa — fa — 1 
formally) , with the exception that we still use the theory for 
/ 2 2 (£i£2 + . . .)o term. 



3.2 Data Vectors 

We intend to analyse the data in the above theoretical frame- 
work using maximum likelihood methods. When construct- 
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Figure 2. This figure illustrates the structure of our covariance 
matrix, which follows from the set up of our data vectors, and the 
strong correlations among different bins of the two- and three- 
point correlation functions. 




Figure 3. The distribution of modes corresponding to the large 
eigenvalues in W in our 22 mock surveys. The mean have been 
subtracted, and each mode has been normalized to unit variance. 
The distribution appears to be consistent with Gaussian. 



ing a data vector from the measurements, we opt for not us- 
ing ratio statistics of the sort Q3 ~ C/£ 2 as it has been done 
in all previous studies, since our measurements include such 
large scales where £ becomes an extremely small number, 
possibly crosses zero. Our ch oice has the addition al advan- 
tage that no ratio bias (e.g., tSzapudi et al.lll999l) appears, 
and our data have direct dependence on as which we will 
exploit. 

We construct the data vector (£, £0) in an appropriate 
scale range, where the theoretical model is expected to be 
a good approximation (c.f. Table |5J. Note that we excise 
redundant entries from the CoCt,^) matrix, i.e. each pair 
of scales appears only once. 



3.3 Covariance Matrix 

In order to perform a \ 2 based maximum likelihood analysis, 
we estimate a covariance matrix from the 22 mock catalogs 
available. These have been extract ed from the Virgo co nsor- 
tium Hubble volume simulation IjColbere et al.ll200(i). For 
detail s of generating catalogues and biasing see ICole et alJ 
(1998). Volume limited subsamples were then created from 
each mock in exactly the same way as we do with the real 
data of 2dFGRS. Thus, for every volume limited subsam- 
ple in Table Q we have 22 measurements of the two- and 
three-point correlation functions. 

We use the usual (biased) estimator for the covariance 
matrix, 



Cij — 



AdfAdj 



(8) 



trix cannot be more than N S i m . In practice we have always 
found this rank to be N s im — 1 in several numerical experi- 
ments. 

To calculate \ 2 horn a singular matrix, we use the sin- 
gular value decomposition (SVD) method to create pseudo 
inverses fc.f. iGaztanaga fe Scoccimarroll200fJ) . We have de- 
composed the C matrix (the tilde denoting estimators is 
omitted from now on) as the multiple of three matrices 

C = UWV T , (9) 

where W is a diagonal matrix, and U, V are orthogonal ma- 
trices; V T means transpose of V. The meaning of this de- 
compo sition is the kern el and image of the linear mapping 
C (c.f. lPress et"allll992fr . Moreover, it is similar to an eigen- 
vector expansion, unique up to degenerate "eigenvalues", 
the elements in W. To calculate x 2 , we need C _1 of our 
singular matrix. Since the inverse is formally VW~ X U T , we 
can replace entries in W 1 corresponding to small eigenval- 
ues with 0; this procedure is usually called constructing a 
pseudo-inverse . 

We have found a marked drop in the eigenvalues beyond 
N S i m — 1. We have checked that the corresponding columns 
of U and V are equal, i.e. they really are "eigenvectors" to 
all intents and purposes. It is meaningful to calculate \ 2 
with the above pseudo inverse, as long as the distribution of 
eigenvectors is Gaussian, as shown in Figure [3 In numerical 
simulations, we found that it is safe to keep about half of the 
eigenvectors, and for all measurements in this paper we used 
the top 10 eigenvalues and their corresponding entries in the 
V matrix as eigenvectors. Our results are robust against this 
choice: have have checked using 8, 9, 10, 11, 12 eigenvectors, 
and the results did not change significantly. 



where Adf = d\ — (dk) is the i-th element of our data vec- 
tor (containing £ and £0 and filtered according the previous 
subsection) from the 7-th simulation; N S i m = 22. 

It can be shown, that a covariance matrix estimated 
by the above equation is necessarily singular if the number 
of simulations is less t han the length of the data vectors 
JSzapudi fc et alJl2005l) . More precisely, the rank of the ma- 



3.4 Maximum Likelihood Analysis 

Using the data vector, and the pseudo inverse of the co- 
variance matrix defined earlier, we can calculate x 2 = 
A T C _1 A, as well as likelihoods p oc exp(— x 2 /2) , where 
A — d ~ d t h(us,b, 62), d t h is the theoretical model with ex- 
plicit dependencies of the parameters. Note that in principle 
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Table 2. Our filtering choices for the data vectors used in the 
analysis; the lower and upper cut offs are displayed in units of 
h Mpc. 3/P stands for 3 parameter maximum likelihood anal- 
ysis based on our phenomenological model, 2/P and 2/T denotes 
2 parameter fits to our phenomenological and theoretical models, 
respectively. 



M b j - 


5 log 10 h 


3/P 


2/P 




2/T 




-18 - 


-17 




11.88 


- 34.97 


18.87 - 


34.97 


-19 - 


-18 


4.04-34.97 


10.18 


- 34.97 


18.87 - 


47.61 


-20 - 


-19 


4.04-34.97 


4.04 - 


34.97 


18.87 - 


64.8 


-21 - 


-20 


5.49-34.97 


11.88 


- 34.97 


18.87 - 


64.8 


-22 - 


-21 




8.73 - 


34.97 


18.87 - 


64.81 



C should also depend on the parameters, but we neglect that 
dependency. This is justified by the final results, which are 
not far from the mock surveys, and simplifies the calcula- 
tion of likelihood enormously. To generalize this method, one 
would need to repeat the simulations and measurements for 
each set of parameters, as well as taking into account the 
determinant in the likelihood; this is clearly unfeasible at 
the moment. 



4 RESULTS 

When applying the above described theoretical framework 
for the interpretation of our clustering measurements, spe- 
cial care needs to be taken about establishing scale ranges 
to be included in the data vectors for maximum likelihood 
analyses. On the one hand, we want to include as much data 
as possible to constrain parameters with the highest preci- 
sion, on the other, if we include data points where the simple 
theoretical model is not a good approximation (for physical 
reasons and/or because of systematic errors), we might bias 
our results. Our considerations are detailed below. 

For the theoretical model, a lower cut of 18h _1 Mpc was 
used, where perturbation theory appears to be very accu- 
rate. In this case we use all scales up to 70h _1 Mpc or the 1/4 
of the characteristic scales of the slice, whichever smaller, to 
avoid severe edge effects. 

For the phenomenological model we use our measure- 
ments from the VLS simulations as theory. Since we neglect 
the errors on these measurements, we use an upper cut of 
35h _1 Mpc, above which the errors are non-negligible. The 
choice of the lower cut is more delicate. There is a com- 
plex interplay between accuracy of the simple bias and red- 
shift distortion models we use with the emergence of dis- 
creteness effects. These finally determine the optimal cut 
in a fairly subtle way. The apparent complexity motivates 
an empirical approach: for the 2-parameter fits we per- 
form maximum likelihood with a series of low cuts between 
4h~ Mpc — 14h~ Mpc for each magnitude limit, and finally 
choose the one with the lowest x 2 . Fortunately we found 
that the two-parameter fits are robust against this choice, 
which nevertheless gives the tightest error bars possible. 

For the three-parameter fit the empirical approach 
would be too expensive and we also want to maximize the 
range as much as possible to resolve the degeneracy between 
as, b. Therefore we use an absolute low cut of 4h~ 1 Mpc, 
below which non-linearities and complexity of the bias are 
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Figure 4. Three parameter likelihood contours for the magni- 
tude slice —21 ~ —20. The wider (blue) curves correspond to 
marginalized distributions, while the narrower (green) curves are 
a cut through the likelihood. The cut through the likelihood curve 
was normalized to the marginalized distribution at the maximum. 
Note that the actual value of the distribution is somewhat arbi- 
trary, as it depends on the grid size. We also show two-dimensional 
marginalized likelihood contours to illustrate the correlations be- 
tween the parameters. 



expected to be strong, complemented with the condition, 
nr 3 47r/3 > 1. This choice gives reasonable control over dis- 
creteness errors which could bias our fit giving a lower cut 
of 5.1h _1 Mpc for the -21 ~ -20 slice. 

In our final choice we took the closest available bin from 
our logarithmic binning scheme to the above values. The 
scale ranges use in the fits are summarized in Table |5] 



4.1 Three-parameter Fits 

We calculated brute force 3 parameter grids with resolu- 
tion 0.005, and 2 parameter grids with resolution 0.001. We 
checked the effects of grid resolution by repeating calcula- 
tions with resolution 0.00025 without any change in the re- 
sults. While erg and b are quite degenerate along the line of 
b x erg = const., as expected, the inclusion of triangles in 
a large range of scales appears to break the degeneracy, at 
least for some of the samples. This is evidenced by the fact 
that a reasonable maximum has developed for the three vol- 
ume limited samples in the mid-range. For those, we present 
our results in Table [3] 

Error bars have been calculated from the marginalized 
curves using 68% thresholds. In the tables, the overall max- 
imum is quoted first with error bars, while the maximum of 
the marginalized likelihood is presented in parentheses. The 
X 2 quoted is normalized to the degrees of freedom, which in 
this case is 10 — 3 = 7. 

The brightest and faintest magnitude limits do not sup- 
port a three-parameter fits. Although they are statistically 
consistent with the three others presented, the large error 
bars due to the degeneracy of erg, b makes them meaning- 
less. Similarly, the theoretical model does not support a sta- 
ble three-parameter fit with either of the magnitude slices. 
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Table 3. Three parameter fits to our two- and three-point clustering measurements. For each parameter, we present the maximum of 
the three-dimensional likelihood with error bars calculated form a 68% of the marginalized distribution. In parentheses we display the 
maximum of the marginalized distribution. 



M bj - 51og 10 h <t 8 b b 2 X 2 

-19 18 l-07tp gi ( °- 80 ) - 81 io'.i2 ( °- 89 ) -°- 06 io.o3 (-0-06) 0.68 

-20 19 0.90lo;°g ( 0.79) OWt^'ll ( 1.01) -0.04l°;^ (-0.04) 3.11 

-21 20 0- 93 ±a20 ( °' 85 ) L04 -o.09 ( 1 - 08 ) -°- 06 ±o.oi (-0-06) 0.97 




1 1.11.21.31.4 %.2 -0.15 -0.1 -0.05 0.05 0.1 



95 1 1.051.11.151.21.251.3 
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Figure 5. Two parameter likelihood contours for the —21 ~ —20 slice. Both the phenomenological model (left) and the theoretical 
model (right) are displayed. 



4.2 Two-parameter Fits 

We also performed two-parameter fits using the bias pa- 
rameters only. These approximately correspond to fitting 
b — > 6(<rg/0.9). Fixing our reference as = 0.9 pins down the 
error bars which would otherwise explode due to the large 
degree of degeneracy, as = 0.9 is consistent with all our 
measurements, and these fits yield extremely precise values 
for the bias parameters. The effective degree of freedom is 
10 - 2 = 8. 

In the case of two-parameter fits, the theoretical model 
also supports stable likelihood surfaces and reasonable error 
bars, with the possible exception of the faintest magnitude 
limit, which we present for completeness only. 



5 SUMMARY AND DISCUSSIONS 

We presented a measurement of the monopole of the three- 
point correlation fu nction in the 2 d FGRS. The new tech - 
nology developed in Szapudi (2004); Pan & Szaoudi (2005) 
enabled the estimation of the latter statistics in a wide range 
of scales from lh _1 Mpc— 70h _1 Mpc. In the three-point func- 
tion, up to 140h _1 Mpc scales enter the measurements and 
the subsequent analysis. In addition, we measured the two- 
point correlation function. 

To interpret these clustering statistics, we developed a 
novel maximum likelihood technique based on joint analysis 
of two- and three-point statistics. We estimated the joint 
correlation matrix from mock 2dF surveys. We seed the top 



10 modes detected by an SVD of the joint correlation ma- 
trix to calculate a generalized x 2 - The distribution of the 
modes is consistent with Gaussian, therefore we maximized 
the Gaussian likelihood with respect the the parameters of 
our theory. Our measurements are robust against keeping 
more or less eigen-modes, as well as against the scale range 
we use for the analysis. In Figure|S|we present the values for 
the bias parameter b as a function of lower scale cut for both 
the phenomenological and theoretical models. Both models 
are fairly robust, as long as the \ 2 1- 

The significance of our three-parameter fit is that it 
yields a highly accurate measurement of as from large scale 
clustering alone. In particular, our estimate is independent 
of the cosmic microwave background (CMB). Yet, our value 
is in excellent agree ment with those der i ved from Wilkinson 
Aniso tropy probe iSpergel et alJ 120031 ; iFosalba fc SzapudU 
|2004|) . as is still one of the most uncertain cosmological pa- 
rameters, and our technique has a great potential to further 
improve the precision of its constraints. The virtually perfect 
agreement of <rg from CMB and large scale three-point level 
clustering is yet another one of the spectacular successes of 
the concordance model. In addition, our result is in excel- 
lent agreement with measurements based on SPS S two-point 
statistics and joint analys is with WMAP llTegmark fc et all 
l2004lPope fe et aljEioI) . 

Our two parameter fits to the bias have extremely small 
error bars, they are likely to be systematics (both in the 
data and in the theory) limited. They can be considered as 
a measurement along the degeneracy line 6crg/0.9 =const., 
and can be directly compared with previous relative bias 
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Table 4. Two parameter fit based on the phenomenological model. Notation is the same as in Table 131 
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2.721 
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Table 5. Two parameter fit based on the theoretical model. Notation is the same as in Table 151 
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Figure 6. We illustrate the degree of robustness of our maximum 
likelihood fit by presenting the best fit bias parameters with error 
bars against the lower limit used in the fit (the upper limit was 
fixed to the value of Table EJ. The (red) squares with error bars 
on the left (scales 4 — 14h — 1 Mpc) show the phenomenological 
model, while on the right (blue) is the theoretical model (scales 
9 — 25h — ^^Mpc). The solid lines show the corresponding x 2 ' s °f 
the fits. In general, as long as \ 2 < 1, the fit is robust. 



mea surements from the 2 dF. Figure[7|presents a comparison 
with lNorberg et alJ i2002l) . and shows excellent agreement. 

We see no statistically significant evidence of scale 
dependency of the bias. Over the full scale range, from 
4h~ Mpc — 70h~ Mpc, a constant bias model gave reason- 
able x 2 ' s - There is no significant trend between the mea- 
surements on intermediate, and large scales using the phe- 
nomenological and theoretical models respectively: the mea- 
sured b's are f ully c onsistent with each other, and with 
iNorberg et alJ J2002T) . The only exception is the theoreti- 
cal model for the faintest subsample, which also has fairly 
degenerate likelihood surfaces. 

A possible check of systematics, although muddled 
by cosmic variance, is to estimate the bias in the NGP 



and SGP s eparately. Our phe nomenological fits corrobo- 
rate that of I Verde etafl teOOST) . who find the SGP slightly 
more biased: our estimate is b = 1.025lo'o48 (1.022), and 
b = 1.0971^22 (1-096) for the NGP and SGP respectively. 
Nevertheless, the two samples are consistent at the 1-1. 5a 
level, supporting the notion that the difference between 
them could be explained with cosmic variance alone. 

The interpretation of our measurements in terms of 
Equations |S|[7| amounts to be some of the most comprehen- 
sive test of our picture of gravitational amplification with 
three point level large scale structure statistics. In partic- 
ular, numerical experiments and theoretical calculations in 
the past tended to focus on a handful of triangular configura- 
tions, with preference for isoceles, and 1 : 2 ratios. We have 
used all possible configurations (although monopole only) 
within our dynamic range and logarithmic binning system, 
and found a good fit to the data. The simplest explanation 
for this is that our basic picture of gravitational amplifica- 
tion is fundamentally correct. In particular, the results lend 
strong support to Gaussian initial conditions, and thus to 
inflation, even though we did not quantify this statement, 
since it was a prior in our model. FigureQshows the remark- 
able success of the theoretical and phenomenological models 
down to lh _1 Mpc. The most natural interpretation of this 
is that bias is relatively simple, and that small scale red- 
shift distortions largely cancel non-linear evolution in red- 
shift space. There is a mild 1.5 — 2a disagreement between 
the theoretical model and the data around 30 — 50h -1 Mpc. 
(c.f. inset of Figure^. While this is not significant according 
to our measured overall \ 2 > it would be interesting to study 
this region with more accurate simulations, and higher order 
theoretical calculations. 

On intermediate scales, we do not detect significant non- 
linear bias. On the other hand, the theoretical model on large 
scales detected non-linear bias at the ~ 2 — 4a level. This 
could mean that either the theory is not accurate enough on 
the largest scales in accordance to the hint provided by the 
inset of Figure Q and/or the largest scales might have some 
yet uncovered systematics, and/or there is significant non- 
linear bias. The latter possibility is somewhat unlikely since 
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Figure 7. Comparison of our bias measuremen ts with the rel- 
ative bias measurements of iNorberg et alJ ^20021) . The solid line 
is 0.85 + 0.15(L/L*), where we normalized L* to be M = -20.5 
magnitude. Red symbols in the middle with large error bars are 
the results of the three-parameter fit. Blue and black symbols 
show the phenomenological and theoretical fits respectively, and 
are shifted to left and right slightly for clarity. Squares and filled 
circles denote the maximum likelihood values in two-dimensions 
and in the marginalized distribution, respectively. 



the intermediate scales do not display significant non-linear 
bias, and it would challenge the well established notion that 
bias should becom e simpler on larger scales. Nevertheless, 
iKavo et alJ (12004) detected a surprizing complexity of non- 
linear bias with the three-point correlation function of the 
Sloan Digital Sky Survey. To decide between the possible ex- 
planations, one would need a highly accurate measurements 
of the three-point function reliable beyond 35h _1 Mpc. While 
this c an and will be done in the future using the Hubble Vol- 
ume llColberg et al . 2000) or similar simulations, the present 
results do not allow distinguishing among the above possi- 
bilities without the risk of over-interpreting the data. 

The technique we presented for constraining cosmologi- 
cal and bias parameters from joint likelihood analysis of two- 
and three-point statistics has enormous potential for further 
high precision cosmological applications. The constraining 
power of the present measurements is limited mainly by the 
theory. For one, only 22 simulations have been used to deter- 
mine the covariance matrix. More accurate covariance ma- 
trix from a much larger number and realistic mocks could 
improve the statistical power of maximum likelihood esti- 
mation based on the same data. In addition, a more realis- 
tic model of bias and redshift distortions of the three point 
statistics, perhaps based on halo models (Ta kada fc Jainl 
120031: iFosalba et al-lliooil . could enable the inclusion of all 
scales measured in the data vector for even tighter con- 
straints. 

The success of the three-parameter fits is remarkable, 
and it is a precursor of potentially even more accurate con- 
straints, and perhaps fits to models with larger number of 
parameters. Perturbation theory can be thou ght of as a 
gener alized bias with an anisotropic kernel (c.f. [Matsubaral 
1995). Since the standard bias model is isotropic, most infor- 
mation on separating bias from gravitational amplification 
should reside in higher order multipoles, such as dipole, and 



quadrupole. Measurements of these in the 2dF will be pre- 
sented elsewhere. In particular, the higher o rder multipole s 
contain information on baryonic oscillations (ISzapudil2 004l. 
which in turn might make it possible to constrain further 
cosmological parameters, such as baryon fraction, and dark 
energy. 

Because of the small number of parameters we used 
so far, a brute force grid technique was feasible. If more 
parameters are fit, our technique lends its elf naturally to 
Mont e Carlo Markov Chain methods (e.g. iLewis fc Bridld 
2002). Along the same lines, a useful and straightforward 
follow up to our investigation is joint likelihood analysis with 
CMB data. This and other generalizations will be presented 
elsewhere. 
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